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In this article we formulate the superperturbation theory for the Anderson impurity model on the real axis. 
The resulting impurity solver allows to evaluate dynamical quantities without numerical analytical contin- 
uation by the maximum entropy method or Pade approximants. This makes the solver well suited to study 
multiplet effects in solids within the dynamical mean field theory. First examples including multi-orbital 
problems are discussed. 
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1 Introduction 

The investigation of strongly correlated materials is one of the main challenges of modern condensed mat- 
ter theory, which has inspired the research activities of many scientists in the last decades. The physical 
properties of such materials are characterized by an interplay between the Coulomb interaction of the 
nearly confined electrons and their kinetic energy. This competition results in a delicate balance of local- 
ization and derealization, which requires accurate non-perturbative approaches. Dynamical mean-field 
theory (DMFT) has become a standard tool for the investigation of strongly correlated materials [1,2] and 
has been applied to models as well as to the investigation of real materials (3]-[5) . DMFT for the first time 
allowed to treat the coherent low-energy excitations and the high-energy excitations as well as their mutual 
feedback on the same footing. 

The main concept of the DMFT is to replace the correlated lattice by a single impurity embedded in a 
self-consistent effective medium. In contrast to classical mean-field approaches, the effective medium in 
DMFT is represented by an energy dependent electronic bath which takes local temporal quantum fluctua- 
tions into account, whereas spatial fluctuations are frozen out. The solution of the DMFT equations in turn 
requires accurate impurity solvers. Quantum Monte Carlo (QMC) algorithms are nowadays widely used 
for this purpose. Continuous-time QMC solvers ||6][7] allow to tackle problems with 5 or even 7 orbitals 
as required for systems with open d or /-shells, respectively. Nevertheless, they suffer from two major 
drawbacks: The fermionic sign problem for a general Coulomb vertex and -more severely- the fact that 
QMC algorithms work in the imaginary time domain require the numerically ill-conditioned analytical 
continuation of stochastic data to the real axis. This makes it difficult to reliably access spectral properties 
and to study e.g. multiplett effects in solids. 

In principal a rigorous approach, which can work on the real- axis, is the exact diagonalization (ED) or 
Lanczos scheme O. Here the continuous bath is discretized by a small collection of bath sites, so that the 
problem can be diagonalized exactly. Due to the exponential growth of the Hilbert space with the number 
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of bath parameters, one is limited to as little as one or two bath sites per d-orbital, so that finite size effects 
may become dominant. In addition, an ambiguity arises in determining the effective bath parameters to 
find an optimal representation of the input hybridization. Alternative approximate approaches comprise 
the iterated perturbation theory (IPT) ||9][T0), the fluctuation-exchange (FLEX) ifTTTl or non-crossing ap- 
proximation (NCA) llT2l and the related hybridization expansion f\3j . These perturbative methods are 
naturally limited to some parameter window and fail outside this region. FLEX and IPT are applicable for 
weak-coupling, while the hybridization expansion works in the strong coupling regime. 
Recently an impurity solver has been developed by combining ED with a diagrammatic approach. The 
key idea is to formulate a perturbation expansion around the ED solution by employing a transformation to 
auxiliary, so-called dual fermions fBl . We refer to the perturbation of a non-trivial (i.e. interacting) albeit 
numerically solvable reference problem as a superperturbation iTHl . It has been shown that the method be- 
comes exact in two opposite limits: For weak coupling and and strong hybridization, the approach becomes 
equivalent to a standard perturbation expansion in the interaction. In the opposite limit of strong interac- 
tion and weak hybridization the formalism resembles the strong coupling expansion |[T3l . This method can 
be sytematically improved by either including more bath sites in the ED (which essentially decreases the 
perturbation) or by including more diagrams. 

In the previous work the solver has been formulated in imaginary time. Analytical continuation was 
found to be significantly more stable than in QMC due to the absence of statistical noise. In this article, we 
formulate the superperturbation theory on the real-axis, making analytical continuation from intermediate 
imaginary time results obsolete. 



2 Superperturbation formalism 

To set the stage, we briefly review the superperturbation formalism in the following. For further reading 
we refer the interested reader to [|T5l . The model under consideration is the Anderson impurity model 
described by the following Hamiltonian: 

H = E 4 7 4At + E e « c « c « + ^[ ct > c i + E (v^lhp + v^% a c p ) . (i) 

&7 ol kot(3 

Here Greek letters are used as a combined index for orbital and spin degrees of freedom, (c^) and b (c) 
are the bath (impurity) creation and annihilation operators, respectively. H mi is the local electron-electron 
interaction, are the transition amplitudes for hopping processes from the bath to an impurity orbital. 
To derive a dual formulation of the problem, we first integrate out the bath degrees of freedom, which leads 
to the conventional action representation: 



S[c*,c] = -}^cZ a [(iw + fi)l-A(o J )Uc^ + S m [c\c}. (2) 

uja(3 

Here S mt is a non-gaussian interaction term and 



fe7 «"-e*7 



is the hybridization function of the full system. The first step in the formulation of the superperturbation is 
to express the action of the model in terms of that of an exactly solvable reference problem and a difference 
term. To this end, we add and subtract a hybridization function A n (cj) corresponding to a discrete bath: 

S[c*,c] = S Kf [c*,c] - ]T 4 a [A N ( W ) - A(lo)Uc^, (3) 
S Kf [c*,c] = - ]T c^Kiuj + - A»] a/3 c w/3 + S int [c*,c}. (4) 

ujo,(3 
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The reference system with a discrete bath with N bath sites shares the interaction part of the original 
problem and can be solved efficiently by exact diagonalization. Note that this step leaves the hybridization 
A n (cj) unspecified. 

The second step is to reformulate the problem in such a way that a perturbative treatment of the dif- 
ference term D(oS) = A N (u) — A(cj) can be performed. Since S ref contains a non-quadratic part Wick's 
theorem is not directly applicable. Therefore, we introduce auxiliary (dual) fermionic degrees of freedom 
using an exact Gaussian identity in the path integral: 

e C* a a a(3 b~^a y sC5 _ \ f X>[f* , J] e -/a 6 «/3//3+/« a «/3C/3+C* dapfp ^ ^ 



where the matrices a and b have the following form: 

b = g-\u)[/^(uj)-A(uj)}- 1 g-\uj) 



CLapb = [A N (o;) - A(uj)] aS = D(u), (6) 



with g(uo) being the exact single particle Green's function of the reference system Eq. ©. After the 
transformation, the resulting action has a mixed representation of c- and /-fermions: 

S[c*,c,f*J] =S re{ [c*,c] + SV*J,c*,c] + J2ti«[9(")D(")9(")]-0U, (7) 
where the coupling between the c and /-fermions is given by 

S C [/*,/,C*, C ] = £affa>)Cutf + £ < a ^H/^- (8) 

uja(3 uja(3 

The original fermionic degrees of freedom can formally be integrated out exactly. To this end, we expand 
in exp(— S c ) in the path integral representation of the partition function. It is convenient to reexpress the 
result in the following form: 



/ 



exp(-S ref [c*,c] - S c [c*,c,r,f})V[c*,c] 

(9) 

= z Ki ex P (- £ f: a g-^)U + V[f*,f}). 

uja(3 



This equation defines the dual potential, which gathers two-particle and higher-order interaction terms. 
The key point is that due to the presence of £ ref on the left-hand- side of ©, integrating out the original 
fermions corresponds to performing the average over the degrees of freedom of the reference system: 
^ref (• • -)ref : = f . . .*D[c* , c]. The resulting dual action has the following form: 

S A [f\f] = -Y.f^[Gl^)]- a lU + V[f\f]. (10) 

As a result of averaging, both the bare dual matrix Green function Gq(uj) = — g(uj)[g(uj)-\-D(uj)~ 1 ]~ 1 g(uj) 
and the dual potential 

V[fi : fi] = ~^Ja{3~/5faf(3f~ff5 + T^lup^S^v fa f(3 fs Jv T • • • • (H) 

contain the correlation functions of the reference system: g(uj) denotes its single-particle Green's function 
and 7^ n ) are the corresponding reducible vertices. The two-particle vertex for example is given by: 

(4) _ -1 -1 r Ref _ 0,Ref i -1 -1 n ~x 

7 a /3 7 5 — 9aa'9 11 '[Xcx'(3' 1 '5' X a '/3'j'5' Mp'pdS'S ' v iZ J 
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Fig. 1 First few diagrams contributing to the dual 
self-energy. Boxes are the reducible vertices of the 
reference system (7^ and 7^), lines denote the dual 
propagator (Go). 

with x ret being the two-particle Green's function of the reference system and ^°' ref its unconnected part. 
These quantities can be calculated straightforwardly on Matsubara frequencies from the Lehmann repre- 
sentation of the single- and two-particle Green's functions |[T5ll . So far, (fTQb is only a reformulation of the 
initial action in Eq. ©. The dual problem can now be treated perturbatively, which essentially corresponds 
to an expansion around the reference problem. For example the self-energy correction stemming from the 
first diagram in Fig. Q]is given by 

Vf a) U = -7$ 75 ( G oW (13) 

After summing up a certain class of diagrams, we obtain a physical result (i.e. solution in terms of the 
physical c-fermion Green function G) by transforming the dual Green's function back to c-fermions using 
the following exact relation [ITU : 

G(u>) = Dicu)- 1 + [g(uj)D(u)]- l G\uj)[D(uj)g(uj)]-\ (14) 

It can be shown that the dual perturbation theory becomes equivalent to conventional perturbation theory 
in the limit of small interaction and strong hybridization. It is instructive to consider its behavior in the 
opposite strong coupling limit. For an expansion around the atomic limit, i.e. A n (cj) = 0, combining 
equations (IT2T) , (TT31) and (ITU) with the lowest order approximation to the dual Green's function, G d w 
^0 + ^o^( a )^o l ea ds to the following expression: 

Gap ~ 9a(3 + g a pf}Tr\gA] - Xap-ysAs-y, (15) 

which recovers the result obtained from an expansion of the imaginary time Green function up to first order 
in the hybridization, as given in Ref. [|T3l . 

The reference system, which is specified through the hybridization function A n (cj), should be chosen 
in an optimal way. Here this is even more crucial than in conventional ED, since for the present approach 
the number of bath sites should be kept at a minimum to increase the efficiency. The choice of A n (cj) 
directly affects the perturbation D(u), which in a certain sense should be minimal. This can be achieved 
by minimizing a predefined distance function ll2lfT6l. 




3 Reformulation on the real axis 

In order to access spectral properties, the previous formulation required analytical continuation of the imag- 
inary time data. This is an ill-posed problem, and different methods have been developed for this purpose. 
The maximum entropy method (MAXENT) ifTTIl was specifically designed for inferring spectral properties 
from statistical data using Bayesian methods. Very fine structures, like multiplets for example, are very 
hard to resolve using this method. In the Pade approach [18] the function of interest is approximated by a 
rational function. Pade can give accurate information on the real axis, provided it is applied to noise-free 
input data. The approximation through a rational function however is unneccessary for the superperturba- 
tion approach and may introduce spurious features. It order to extract unbiased real axis information, it is 
desirable to formulate this approach directly on the real axis. 
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In ED, real axis information is readily obtained by performing the substitution iuo —> E + iS in the 
Lehmann representation, where S is a small broadening parameter. For the single-particle Green's function 
this leads to: 

n,m 

From latter expression the density of states (DOS) is obtained as DOS(E) = — l/irImg(E + In 
the superperturbation, analytical continuation is more involved, as one has to take care of the individual 
diagrams. In the following we discuss the analytical continuation of the first diagram (diagram a.) of Fig. 
[]])• This diagram has been found to yield by far largest correction to the initial solution. Inserting (Tf2l) into 
JT3l L we obtain 

= -5^/3Tr[A d 5 ] + A d a/3 - g^Xa'P'YS'gJp^y, (17) 

which is very similar to expression (H5T l where now A d := [g + (A N — A) -1 ] -1 plays the role of A. 

While g and A N are readily accessed on the real axis , we can make use of a result of Ref. fT3l which 
allows to calculate the product of the two-particle Green function with A d directly within ED. 

ioo n / 7<5 ijkl 

./ * * x [ 7Zjs(Ej,Ek) ! 1Z 7 s(Ei,Ek) ! Q 1 8(iu,E i ,Ek) 1 



Eji(iuj - Eji) Eij(iuj - En) (iu - E u )(iu - Eji)- 
lZ 1 5(E k , Ej) TZ^sjEk.Ei) Q 1 s(—iuj,E k ,Ei) - 

1 l E jt (iuj - E^) + - £ Zi ) (iw - £ m )(«j - 

IZrysjEj^ Ej) TZ^sjEj, Ej) Q 1 s(—iuj, E^ Ej) 



+ (c a CsC*Cp)ijkl [ 
+ (c*CsC a Cp)ij k i [ 
+ {c5C*C a Cp)ijkl [ 



E ki (iuj - Ei k ) E ik (iuj - En) (iuj - Eu)(iuj - E tk )- 
1Z 7 s(Ejj Ek) lZ 7 s(Ej : Ei) Q 7 s(iu, Ej, Ei) 



(18) 



-E ki (iuj - Ei k ) E ik (iuj - En) (iu - Eu)(iu - E tk ) j 

7 p/ (lUJ - E kj ){lUJ - En) 

[n 75 (Ek,Ei)-n 7 5(Ej,E z ) + Q l5 (iu,Ej,Ei) - Q l5 (-iu, E k , Ei)] 

7 P {iu ~ E kj )[iuj - En) 
\jZ 1 s{Ei,E k ) — lZ 1 s{E i ,Ej) + Qrys(iuj,Ei,E k ) — Q 7 s(—iv,Ei,Ej)\, 

with the following definitions for the matrix elements: 

(O a OpO y O s )ijki = (i\O a \j)(j\Op\k)(k\0 7 \l)(l\O s \i). (19) 
The functions 1Z and Q have the following definitions: 

^<«,^ !(.-«+.-«<) ,20, 



r_| e -^ A d M for Ei = E d 

Qcpi^EuEj) := < —BE \ 1 (2D 
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Fig. 2 top: Standard procedure of 
converting a Matsubara sum into a con- 
tour integral. In a first step the sum is 
replaced by a contour integral over the 
function itself times the Fermi func- 
tion, which has poles at Matsubara fre- 
quencies. In the second and third step 
the contour is deformed to two line in- 
tegrals along the real axis. The illus- 
trated contour deformation is applied 
in the calculation of 1Z. hot. : Construc- 
tion of the contour integral to calculate 
Q. In comparison to the upper case one 
has to take into account an additional 
pole at z = iuo + En- 



The trace over A D and g can also be expressed in terms of 1Z\ 
Tr[ 5 A d ] = (i\c S \j){j\c\\i)K l5 {E^E 3 ). 



(22) 



Hence the analytic continuation of Eq dT71) reduces to the continuation of the function Q(iuo, E iy Ej). Since 
1Z(Ei,Ej) does not depend on iuo the function can be in principle calculated using definition (l20l) . In order 
to treat both function evaluations on the same footing, we discuss in the following how 1Z and Q can be 
calculated via an integral along the real axis. 

The Matsubara sums in Q and 7Z are rewritten as a sum over residues of the function itself times the Fermi 
function, which has poles of first order at the Matsubara frequencies: 



±YF(ico) = -^- 



2iti J c e 



F(E) 



dE. 



(23) 



Afterwards the integral contour is deformed according to figure This leads to the following definitions 
of both expressions as integrals along the real axis: 



n{E u E 2 ) 



X x +X 2 



2ni 



Q(iu,E 1 ,E 2 ) =(X 1 -X 2 ) 



f 

J — c 



A i (z+)f(z+) 



J — ( 



A d (z-)f(z-) 



' A d (^+)/(z+) 



E 



12 



dz, 



oo 



y + 



f 

J — c 



A«(z-)f(z-) 



E 



dz 



12 



iuo — E\2 
1 



dz 



1 



oPEi 



■A d (iuo + E 12 ) 



(24) 



(25) 



where z ± = z±ie, withe < ir/(3 being the offset of the contour from the real axis and Xi = ex.p(—f3Ei)/Z. 
The last term in the expression for Q is due to a residue of Q at z = iuo + £'12- Now the analytical contin- 
uation can be performed by replacing iuo by E + iS: 



Q(z,E 1 ,E 2 ) = (X 1 -X 2 ) 



1 



30 A d (z+)/(z+) 

**(*->/(*-> dz 



o/3E! 



A a (z + E 12 ) 



(26) 



© 2007 WILEY- VCH Verlag GmbH & Co. KGaA, Weinheim 



www.ann-phys.org 



Ann. Phys. (Berlin) 0, No. (0) / www.ann-phys.org 



7 




Fig. 3 SPERT calculations on the real axis 
(CORA), in comparison with Pade and the 
result of the reference system (ED). The 
upper part shows an example for f3 = 
10 and U — 3, in units of the half- 
bandwidth (W/2). Here the agreement be- 
tween CORA and Pade is good, but the 
CORA shows more structure on the real 
axis. The lower part shows the same system 
with an applied magnet field: B = 0.05. 



Fig. 4 Example of a CORA calculation 
for a multi-orbital case. A three orbital im- 
purity with U = 1.5, U' = 0.7, J = 0.4, 
H — 3.02 and ft = 15 (in units of the 
quarter-bandwidth) is embedded in a bath 
with a constant density of states in the en- 
ergy window \E\ < W/2. The coupling 
to the bath is V = 0.1. The hybridiza- 
tion was approximated by a single bath site. 
The CORA data set exhibits small features 
in the density of states which are absent in 
the Pade results. 



with z = E + iS, where 5 is the usual broadening parameter, which is restricted to values S > e. This 
completes the analytic continuation. 

A few remarks are in place. In contrast to Ref. [fT3l we do not perform the limit e —> because of the 
particular structure of the "dual" hybridization, which has poles on the real axis. The integrals are evalu- 
ated with a small offset e. For not too low temperatures a simple quadrature rule is sufficient. For lower 
temperatures, A d develops a sharp peak. We therefore employ an adaptive Gauss-Kronrod algorithm taken 
from the GNU scientific library ||T9l . As a check for numerical accuracy one may verify that the value of 
the integrals is independent of e. 

Figure[3]shows some illustrative results for the AIM with hybridization corresponding to a semielliptical 
DOS of bandwidth W. In the upper left plot the calculation on the real axis (CORA) is compared to an 
analytic continuation using Pade and the initial solution of the reference system (labeled 'ED'). The ED 
curve has a clear splitting at the Fermi level, whereas the CORA curve exhibits a Kondo peak. The CORA is 
in a good agreement with Pade. The lower left plot shows an additional example with an applied magnetic 
field. Here the splitting of the peaks is clearly visible and the CORA is again in good agreement with Pade. 
For completeness the data on Matsubara frequencies has been added on the right. 

In Figure [4] we present an example for a multi-orbital problem with a rotational invariant Coulomb vertex. 
A three-orbital impurity has been embedded in a bath, which corresponds to a flat density of states in the 
energy window \E\ < W/2. The coupling to the bath was moderate and has been approximated by a 
single bath site, which was equally connected to all impurity orbitals. In comparison to the solution of 
the reference system a clear shift of nearly all peaks in the CORA is visible. The Pade solution is in good 
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overall agreement with the CORA data, but fails to reproduce some small features. Indeed one may expect 
that the CORA results are more accurate than the analytical continuation via Pade approximants if the 
continued function has rich structure, since the approximation through a rational function becomes less 
accurate. For complicated multiorbital systems we expect differences to be more pronounced. 

4 Conclusion 

In the present work we have presented a multiorbital impurity solver based on the superperturbation for 
the Anderson impurity model (AIM). It allows to compute dynamical quantities directly on the real axis, 
which has the advantage that no analytic continuation using approximate methods like MAXENT or Pade 
is necessary. Considering three examples including a multi-orbital impurity problem, we compared the 
results of the new implementation to results obtained using Pade. We find overall good agreement. The 
direct calculation on the real axis however can resolve finer structures, which will be useful for the study 
of multiplett effects in solids. This work has been supported by the Cluster of Excellence Nanospintronics 
(LExI Hamburg) and DFG Grant(436113/938/0-R). 
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